Canonical Bose gas simulations with stochastic gauges 
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A technique to simulate the grand canonical ensembles of interacting Bose gases is presented. Re- 
sults are generated for many temperatures by averaging over energy-weighted stochastic paths, each 
corresponding to a solution of coupled Gross-Pitaevskii equations with phase-noise. The stochastic 
gauge method used relies on an off-diagonal coherent-state expansion, thus taking into account all 
quantum correlations. As an example, the second order spatial correlation function and momentum 
distribution for an interacting ID Bose gas are calculated. 
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Calculating the observables corresponding to a grand 
canonical density matrix is a highly nontrivial problem in 
quantum many-body theory, due to the large dimension- 
ality of the corresponding Hilbert space. In this Letter, 
we show that for Bose gases, this problem is soluble using 
a coherent-state stochastic gauge method, which gener- 
ates equations similar to the widely used Gross-Pitaevskii 
equation Q], with additional quantum noise terms. The 
issue is of much topical interest, as it relates not just to 
recent experiment and theory on trapped dilute gas Bose- 
Einstein condensates (BEC) [2|, but also to fundamental 
questions like the fermion-boson duality problem Q in 
quasi one-dimensional (IP) Bose gas systems 

A major advantage of the method presented is that 
the number of variables is linear in the number of spatial 
lattice points, in any number of dimensions. This scaling 
suggests it is competitive in efficiency with Monte Carlo 
techniques, while being able to in principle calculate ar- 
bitrary observables. Additionally, it is applicable to both 
imaginary time evolution equations describing thermal 
equilibrium and to quantum dynamical calculations in 
real time. As a result, the same technique can be used to 
calculate state preparation and subsequent time evolu- 
tion. In this Letter we calculate equilibrium momentum 
distributions and spatial density correlation functions for 
a finite-temperature IP Bose gas, as well as verifying the 
technique by comparisons with known exact results. 

The exactly solvable ID Bose gas model with repul- 
sive delta function interactions has been the 
subject of renewed theoretical interest recently, owing to 
experimental realization of the one-dimensional regime 
in trapped alkali gases 0|. The higher-order correlation 
functions are responsible for the rates of inelastic colli- 
sional processes, and are of particular importance for the 
studies of interference and coherence properties of atom 
lasers operating in one-dimensional waveguide environ- 
ments. Local higher-order correlations at zero and finite 
temperatures have recently been calculated 0,0, utiliz- 
ing the previously known exact analytic solutions 

Here we employ the stochastic gauge method to obtain 
results for the non-local second-order correlation func- 



tion. In particular, we investigate the signatures of the 
cross-over from the weak-coupling Gross-Pitaevskii (GP) 
regime to the strong-coupling Tonks- Girardeau (TG) 
regime where the system behaves like a free Fermi gas. 
The results indicate that these correlations are a stronger 
indicator of quantum effects than the recently observed 
momentum distributions (lfij . 

The stochastic gauge method is a generalization of the 
positive P -representation approach, and is described in 
detail in Ref. The original P-function expansion 

in quantum optics using coherent states^], was due to 
Glauber and Sudarshan I3|. This does not give a pos- 
itive distribution for all density matrices and was later 
generalized to non-diagonal P-functions fliL llfij. which 
always exist and are positive. 

This positive P-representation has been successfully 
applied to quantum dynamical problems in quantum op- 
tics 0| and many-body theory 0|. If there is no damp- 
ing, the method has stability problems which are 
overcome by introducing stochastic gauges Other 
known methods include the quantum Monte-Carlo meth- 
ods and stochastic wave-function methods [2(1 
which have been applied to BECs. However, these cannot 
be used for grand canonical ensembles, nor for subsequent 
quantum dynamical calculations if there are losses. 

We start by considering the interacting Bose gas model 
with repulsive delta-function interaction between the par- 
ticles. In second quantization, the model Hamiltonian is 



H = 



U 



+V(x)^(x)^{x) + — tft(x) 2 tf(a;) 2 



dx. 



(1) 



Here ^(x) is the field operator, m is the mass, U > 
is the coupling constant for the inter-particle interaction, 
and V(x) is an external potential. 

Restricting ourselves to the ID case, with the under- 
standing that U corresponds to the IP; coupling strength 

E3| ; let us reduce Eq. QJ to a Bose-Hubbard lattice 
Hamiltonian 1231 which contains all the essential features. 
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For a lattice with M sites and spacing Ax, defining lat- 
tice annihilation operators a — {dj} = {V Ai$(ij)}, one 
obtains (summing over repeated indices): 



X 



-2u>ija\a,j 



(2) 



and Uj = djdj 



In this Hamiltonian, the cjy represent the dimensionless 
oscillator frequencies together with a linear coupling to 
other sites, while x = U/2Ax is the rescaled coupling. 
The boson commutators are [fij,aj 
is the boson number operator for the j-th site. 

To extend the gauge P- representation equations ^l] 
to imaginary time or temperature, consider that when 
[H, N] = , the un-normalized density matrix of the 
grand canonical ensemble is given by 



Pu 



-(H~[j,N)/k B T 



-Kt 



(3) 



Here fi(T) is the temperature-dependent chemical poten- 
tial, N is the particle number operator, and we introduce 
dimensionless 'Kamiltonian' K{p,aj,c^) = (H — pN)/x 
and r = x/(ksT) corresponding to inverse temperature. 
In the Schrodinger picture, Eq. <EJ| leads to the "imagi- 
nary time" master-like equation 



dp u 
dr 



iPu 



(4) 



where p e = 9[r/x(r)]/9(2xr) and [A, B]+ stands for anti- 
commutator. The usefulness of this equation relies on the 
fact that at initial "time" r = (i.e. T — > oo) the state 
is p u (0) — exp(— AiV), where A = — lim r ^o[ r A i ( T )]/x ) 
and the initial number of particles per site is no = 
[exp(A) — This initial state can then be evolved in r 
to obtain grand canonical ensembles at all other (lower) 
temperatures T. 

We expand the density matrix with a distribution G 
and an additional quantum amplitude fl: 



Pu 



4M+2- 



(5) 



where A = D, \\a) (/3*|| e~ a/3 , a and f3 are complex M- 
dimensional vectors of Bargmann coherent state ampli- 
tudes with elements at each lattice point, and the no- 
tation ~a is shorthand for all the variables ($l,a,(3) = 
(ao, ai, . . . , o>2m)- The initial gauge-distribution is: 

G (a) oc exp [- \a\ 2 /n ] S 2M (a - (3*)S 2 (n - 1). (6) 

To determine the effects of the time evolution <@J we 
use standard differential identities 0] to replace the an- 
nihilation and creation operators acting on the projector 
A . Using these operator identities the operator equation 
(@J can be transformed to: 



4M+2- 



(7) 



We can define the resulting differential operator L in two 
parts as C = C K + C 9 , by using the operator identities to 
obtain C K = -K(t) + ffj=i [Afdj - a]d]/2] - and we 
will define the "gauge" term C 9 next. Here, Af — (fi e — 

nAoti —LUijaj, U)i+M,j+M = ^ji,<^i+M,j — <^i,j+M = 0, 

and dj — d/daj. The effective c-number 'Kamiltonian' 
is K(t) = K(2xPc{t) ,a,/3), while the effective complex 
boson number is rij = n'j + in" = n.j + M = otjPj- 

Now let us introduce a suitable stochastic gauge oper- 
ator C 9 . The do identity 0] tldoA = A allows the use 
of arbitrary stochastic gauge functions. This is similar 
to an electrodynamic gauge, in that it defines an infinite 
class of physically equivalent Fokker-Planck and hence 
stochastic differential equations. To show how this is 
achieved, define 2M arbitrary complex gauge functions 
g = {gi(~a , a *)} to give a gauge-dependent differential 
operator C 9 which satisfies C 9 A = 0: 



C 9 = 



-k(t) + ls 2 nd { 



1M 



(ttd - 1) 



(8) 

With a suitable gauge choice that eliminates bound- 
ary terms, we can now transform the equation Q into a 
positive-definite (ill fl4| Fokker-Planck equation via par- 
tial integration, and then into Ito stochastic differential 
equations using standard methods 



^± = A j + ia j Q(t) 



2M 
i=l 



(9) 



The total complex drift vector including the quantum 
amplitude Q, is now A = (— QK(t), A%, . . . , A2m)> where 



for j > 0, 



A 



K 



igjctj. Here Co = 0, and 



Cj(t) are 2M independent Gaussian noises such that 
(Ci(rMr'))=S ij S(r-r'). 

With no gauge (gi = 0), one finds that for this sys- 
tem the equations are unstable for n\ < because they 
contain terms of the form d, = —a 2 /3i. These have been 
shown to lead to systematic errors due to singular tra- 
jectories which cause boundary terms 0|. We choose 
gauge functions of the form gj — i(n'- — \nj\), which re- 
move the instabilities, and also stabilize the phase of rij 
so that (3j ~ a*. This then gives the stable equations 



da 



2M 



= [ji e - | rij | - in"} ^ - ^wya^ + jaiCi(r), 



dr 

dtt 
~dr 



2M 



n. 



(10) 



There is an intuitive physical interpretation. Since — 
a* in the initial thermal state, each amplitude initially 
obeys a Gross-Pitaevskii equation in imaginary time, 
with quantum phase-noise due to the interactions. This 
causes non-classical statistics with a ^ f3* to emerge as 
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the temperature is lowered. Along each path an addi- 
tional ensemble weight fl is accumulated, which is pro- 
portional to the effective Gibbs factor e~ J K ( T ) dT for the 
path, together with gauge-dependent terms. This gives 
a hybrid between a path integral and a purely stochastic 
type of simulation. 

Observables, which are always expressible as sums of 
terms of the form aJ n a m , have expectation values given 
by the stochastic averages 



(n + n* 



stoch 



and can all be in principle calculated. 

As a demonstration, we present the results of a uni- 
form (untrapped) ID Bose gas calculation in a regime 
which lies in the transition region between the weakly- 
and strongly-interacting gas, where neither perturbation 
theories nor the TG Fermi gas approximations work well. 

The behavior of a uniform Bose gas described by the 
model (QJ in ID depends on two parameters T le \ = T/Td 
and 7 = mll/ph 2 . For an ideal gas with linear (ID) par- 
ticle number density p, the transition from Boltzmann 
to Bose statistics occurs around the quantum degener- 
acy temperature T c i = 2irti 2 p 2 /mkB- This corresponds 
to the temperature when the average spacing between 
particles l/p equals the thermal de Broglie wavelength 
A t h = \j2irh 2 /mkBT . 

The second parameter of interest is the coupling 
strength 7. As 7 increases from zero, one moves from 
the weakly-interacting GP regime, where the gas is well 
described by Bogoliubov theory, to a strongly interact- 
ing TG regime, where it undergoes "fermionization" 0,0- 
For example, at zero temperature the ground state en- 
ergy is described Q to within ~ 95% accuracy by Bo- 
goliubov theory for 7 < 0.8, and by an ideal Fermi gas 
for 7 > 80. Here we simulate for an intermediate regime 
with 7 = T re ] = 10, which is more likely to be accessible 
experimentally. 

A point to note is that the chemical potential of the 
diffusive reservoir in contact with the system /i(r), can 
be chosen at will. Its form must fulfill two conditions: 

(1) The desired values of 7(T re i) are simulated. Since \ 
is a constant during the simulation, /i(r) must be tailored 
to give the desired densities p(r), and hence 7. 

(2) A finite lattice leads to a maximum momentum 
cutoff, so one must take care that this does not influ- 
ence observables — i.e. that simulated momenta do not 
couple to momenta beyond the cutoff. Fortunately at 
high enough temperatures T re i ^> 7, kinetic processes 
dominate over interparticle interactions and the momen- 
tum modes decouple. Hence, it suffices to start the 
simulation at To in an ideal Bose gas state at a high- 
enough temperature so that T re i(ro) ^> 7 (to). The ini- 
tial state is then still given by a form like (HI, but the 
initial occupations in momentum space no(k) are pro- 
portional to the ideal gas Bose distribution no(k) oc 




FIG. 1: ID Bose gas in an intermediate regime: T re i = 10 
and 7 = 10. Left panel shows the momentum distribution, 
while the right panel is the second-order spatial correlation 
function. Here, the distance x is in units of the he aling length 
f, which can be expressed as £ = A t h^/2T re i/27. Solid lines 
are the results of numerical simulations with 738200 stochastic 
trajectories. Dashed lines are the results for 7 = ideal Bose 
gas (ID), 7 — » 00 Tonks-Girardeau (TG) gas, and T re i — > 00 
Boltzmann (BL) gas, shown for comparison. 



{exp[(— h 2 k 2 /2m + p,)r /x] — I} -1 - In the example con- 
sidered here T le ](T )/-f(r ) ~ 2.01 x 10 3 and this ratio is 
monotonically decreasing with r. 

To avoid wasteful calculation of empty modes in the 
simulation, the width of the momentum distribution 
of PyTiei should be kept approximately constant 
throughout. Hence, we want p(p(r)) oc 1/ y/T ie \(r). A 
convenient form of p,(r > t ) gives a smooth interpo- 
lation between an initial fugacity zq = cxp(tqp,(tq)/x), 
and a final fugacity Zf. We choose the function: 2z(t) — 

Z + Zf + (zj - Zq) COS (tt[t - Tf]/[Tf - T ]) . 

To obtain the desired 7 = T re i = 10 at the target tem- 
perature Tf — x/^BTf [criterion (1)], one must choose 
two of the three numerical parameters x, Tf, and Zf, while 
the third becomes a scaling parameter which determines 
the choice of units. To also fulfill criterion (2) one must 
choose such values of zq and To that T Te \(To, zq, x) 3> 
7(to,zo,x)- We chose the values zq — Z//1000 and 

T = T//4. 

Figure compares two important observables for the 
7 = Trei = 10 homogeneous system, with those of ideal 
Bose (7 = 0) and TG (7 — > 00) gases under the same 
temperature and chemical potential. A comparison is 
also made to the Boltzmann distribution which occurs 
at high temperatures T re \ — > 00. The particle number 
density in momentum space (classically of fc-space width 
kth = Ath/V^r) is intermediate between the ideal gas 
and the TG regime. At high momenta, all the tails 
are seen to be Boltzmann-like, whereas for low momenta 
Bose enhancement occurs above the Boltzmann level, but 
to a much smaller degree than in the ideal gas. 

The second-order spatial correlation function 



5 (2) (*) = 



(Wt(Q)^t( a; )^( a ;)^(0)) 

(*t(0)*(0))(#t( x )«f(V)) 



(12) 
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FIG. 2: Comparison of numerical calculation (solid lines) to 
exact results [3 (asterisks), as well as to Ideal (ID) and Tonks- 
Girardeau (TG) gases at the same temperature T and chem- 
ical potential. Left panel: particle number density. Right 
panel: energy per particle £ . 



quantifies the spatial clumping of particles. A spatially 
uncorrelated field has g^ 2 \x > 0) = 1. Some features 
that can be seen from Fig. include: (i) A peak arises 
at x ~ 0.5 £, corresponding to the most likely separation 
between particles. This effect is not present in an ideal 
gas. (it) Quite strong anti-bunching (g^ 2 '(0) < 1) of the 
bosons is seen: 3 (2) (0) = 0.72 ± 0.01. (in) While the 
momentum distribution p(k) is quite similar to that of 
a Boltzmann gas, the behavior of the pair correlation 
function g^ 2 \x) is widely different. 

As an added test, we compared these calculations with 
exact results for boson density and energy derived by 
Yang and Yang 0. This comparison is shown in Fig. 
El and one sees that the agreement is excellent. Similar 
good agreement is found with the value of g^ (0) calcu- 
lated PI from the exact Yang- Yang theory. 

In summary, we have applied the gauge P- 
representation technique to derive a set of stochastic 
equations whose averages allow the calculation of all ob- 
servables for the grand canonical ensemble of an inter- 
acting Bose gas modeled by Eq. QJ. Results are ob- 
tained for a range of temperatures in one simulation. In 
this, an ensemble of trajectories representing a quantum 
many-body system is cooled from an initially high tem- 
perature. The technique is readily applicable to more 
dimensions and external potential by a straightforward 
generalization of the derivation outlined here. Further 
improvements are possible through optimizing the ba- 
sis (e.g. with generalized Gaussian states), and by us- 
ing other gauges and algorithms such as Metropolis sam- 
pling. 

We have used this method to calculate the second- 
order spatial correlation function and momentum distri- 
bution for an interacting ID Bose gas, in a regime where 
both strong and weak coupling approximations are in- 
valid, and no exact results are currently available. Our 
results show that the pair correlation function is a sensi- 
tive indicator of the Bose-Fermi crossover, even at tem- 
peratures above quantum degeneracy. 
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